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In this paper we analyze the performance of the Quantum Adiabatic Evolution algorithm 
on a variant of Satisfiability problem for an ensemble of random graphs parametrized by 
the ratio of clauses to variables, 7 = M/N. We introduce a set of macroscopic parameters 
(landscapes) and put forward an ansatz of universality for random bit flips. We then for- 
mulate the problem of finding the smallest eigenvalue and the excitation gap as a statistical 
mechanics problem. We use the so-called annealing approximation with a refinement that 
a finite set of macroscopic variables {versus only energy) is used, and are able to show the 
existence of a dynamic threshold 7 = 7^ starting with some value of K - the number of 
variables in each clause. Beyond dynamic threshold, the algorithm should take exponen- 
tially long time to find a solution. We compare the results for extended and simplified sets 
of landscapes and provide numerical evidence in support of our universality ansatz. We 
have been able to map the ensemble of random graphs onto another ensemble with fluctu- 
ations significantly reduced. This enabled us to obtain tight upper bounds on satisfiability 
transition and to recompute the dynamical transition using the extended set of landscapes. 

PACS numbers: 03.67.Lx, 89.70.+C, 05.20.-y 



I. INTRODUCTION 

An important open question in the field of quantum computing is whether it is possible to 
develop quantum algorithms capable of efficiently solving combinatorial optimization problems 
(COP). In the simplest case the task in a COP is to minimize the energy function E a with the 
domain given by the set of all possible assignments of A" binary variables, <r = joy, . . . , a N }, 
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(Tj = ±1. In quantum computation this cost function corresponds to a Hamiltonian Hp 

Hp = Y^E<r\a)(<T\ (1) 

|cr) = \a x )i (8) |<7 2 ) 2 ® • • • ® |o"at)at, 

where the summation is over the 2 N states |cr) forming the computational basis of a quantum 
computer with iV qubits. State \<Jj)j of the j-th qubit is an eigenstate of the Pauli matrix a z with 
eigenvalue <jj. It is clear from the above that the ground state of Hp encodes the solution to the 
COP with cost function E„. In what follows we shall use two equivalent notations for binary 
variables: Ising spins cr,- = ±1 as well as bits Zj = (1 — (Jj)/2 = 0, 1. 

Recently Farhi and coworkers proposed a new family of quantum algorithms for combinato- 



rial optimization that is based on the properties of quantum adiabatic evolution QJJ, |2D - Numerical 
simulations were performed for the study of its performance for satisfiability problems [3]. Im- 
plementation of these algorithms on a quantum computing device is feasible for COPs where the 
energy function E„ possesses a locality property, in a sense that it is given by the sum of terms 
each involving only a relatively small number of bits, that does not scale with iV III 0, If]]. An ex- 
ample of a problem that can have this property is Satisfiability that deals with iV binary variables, 
submitted to M constraints, assuming that each constraint involves 0(1) bits. The task is to find a 
bit assignment that satisfies all the constraints. 

Sa, isfl a bilily is a basic prohien, ,„ the so-caHed NP-co mp ie,e class Q. This class co„,a,„s 
hundreds of the most common computationally hard problems encountered in practice, such as 
constraint satisfaction and graph coloring. NP-complete problems are characterized in the worst 
case by exponential scaling of the run time or memory requirement with the problem size N. A 
special property of the class is that any NP-complete problem can be converted into any other 
NP-complete problem in polynomial time on a classical computer. Therefore, it is sufficient to 
find a deterministic algorithm that can be guaranteed to solve all instances of just one of the NP- 
complete problems within a polynomial time bound. It is widely believed, however, that such an 
algorithm does not exist on a classical computer. Whether it exists on a quantum computer is one 
of the central open questions. 

Running of the quantum adiabatic evolution algorithms (QAA) for several NP-complete prob- 
lems has been simulated on a classical computer using a large number of randomly generated 
problem instances that are believed to be computationally hard for classical algorithms J^lslQ]. 
Results of these numerical simulations for relatively small size of the problem instances ( N < 
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25) suggest a quadratic scaling law of the run time of the QAA with N. 

A particularly simple version of Satisfiability is the NP-complete Exact Cover problem that 
was used in [2] to study the performance of QAA. In this problem each constraint is a clause that 
involves a subset of K = 3 binary variables. A given constraint is satisfied if exactly one of its bits 
equals 1 and the rest of the bits equal 0. In the optimization version of this problem one minimizes 
the energy function that is equal to the number of constraints violated by a given bit-assignment 
cr. A generalization of this problem to an arbitrary number K can be called positive 1-in-K SAT 

Q. 

In practice algorithms for NP-complete problems are characterized by a wide range of running 
times, from linear to exponential, depending on the choice of certain control parameters of the 
problem (e.g., in Satisfiability it is the ratio of the number of constraints to the number of vari- 
ables, M/N). Therefore, a practically important alternative to the worst case complexity analysis 
is study of a typical-case behavior of optimization algorithms on ensembles of randomly generated 
problem instances chosen from a given probability distribution. For example, in the case of pos- 
itive 1-in-K SAT one can define a uniform ensemble of random problem instances. An instance 
X consists of M statistically independent clauses, each corresponding to a if-tuple of distinct 
bit-indices uniformly sampled from the interval (1, N) with probability 1/ (^) . 

In the case of an exponential scaling low for the algorithm's running times t a it is convenient to 
analyze the distribution of a normalized logarithmic quantity \ogt a /N. This distribution becomes 
increasingly narrow in the limit of large N where the mean value (log t a ) /N well characterizes the 
typical case exponential complexity of an algorithm. For Satisfiability problem the dependence of 
the asymptotic quantity 

77 = lim {\ogt a )/N (2) 

on the clause-to-variable ratio 7 = M/N has the qualitative form shown in Figd At some critical 
value J — jd algorithmic complexity undergoes the dynamical transition from polynomial to 
exponential scaling law. This transition has been studied recently for the case of a variant of 
the classical random-walk algorithm for the Satisfiability problem [9]. Function 7/(7) is non- 
monotonic in_7 and reaches its maximum at a certain point 7 C > jj. It was discovered some time 
ago II2II that 7 C is a critical value for the so called satisfiability phase transition: if 7 < j c , 

a randomly drawn instance is satisfiable with high probability, i.e., there exists at least one bit 
assignment cr that satisfies all the constraints (E a = 0). For 7 > 7 C instances are almost never 
satisfiable. In the asymptotic limit N — > 00 the proportion of satisfiable instances drops from 1 to 
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FIG. 1: Solid line shows the qualitative plot of the normalized quantity rj/r] maiX vs M/N (r/ max is a maxi- 
mum value of rj). Dashed line shows the proportion of satisfiable instances vs M/N. 

infinitely steeply at 7 = 7 C as shown in Fig. [I] 

The value of 7^ (unlike 7 C ) depends on both the problem at hand and the optimization algo- 
rithm. Recent years have seen a rising interest in study of dynamic threshold phenomena for local 
search algorithms J9 , [13I] . That effort is in its initial stage and simple approximations (in spirit of 
annealing approximation) were employed to estimate the location of threshold. Comparison of the 
dynamical thresholds 7^ for different algorithms provides an important relative measure of their 
typical-case performance in a given problem. 

This paper is organized as follows. In section|n|we introduce the Quantum Adiabatic Evolution 
Algorithm and explain how the complexity of the algorithm depends on the spectrum of the Hamil- 
tonian. In sectionHTIlwe formulate quasiclassical approximation used to study the complexity and 
introduce the notion of landscapes. In section[[V]we introduce positive 7\T-NAE SAT and positive 
l-'m-K SAT - the NP-complete problems, which we use as a test bed for our method. In section 
IVl we provide detailed computation of entropy and landscapes within annealing approximation. 
We discuss the universality of landscape probability distributions in section|Vl| Sections fVIII and 
IVIIII are devoted to improving the annealing bound. A subgraph responsible for the hardest part 
of the problem (a core) is identified and results are rederived for the subgraph. In all cases we 
are concerned in finding the dynamic threshold - the critical ratio of clauses to variables above 
which the algorithm is expected to take exponentially long time to find a solution. We discuss 
our results as well as possible ramifications and extensions of our work in Conclusion (section 
[Txli. Appendix lAl gives a sketchy proof of NP-completeness of the problems we considered, and 
appendix |H1 discusses an incremental improvement over annealing approximation possible within 
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our formalism. 



II. QUANTUM ADIABATIC EVOLUTION ALGORITHM 

Consider the time-dependent Hamiltonian H(t) = H(t/T) 

H(t) = (1-t)H b + tH p , (3) 

where r = t/T E (0, 1) is dimensionless "time", Hp is the "problem" Hamiltonian ([TJ and Hp 
is a "driver" Hamiltonian, that is designed to cause transitions between the eigenstates of Hp. 
Using dimensionless time and setting h = 1 the quantum state evolution obeys the equation, 
iTd\^(r))/dr = H(t)\^(t)). At the initial moment the quantum state |^(0)) is prepared to be 
the ground state of H(0) = Hp- In the simplest case 

N 

Kb = -J2< |v&(0))=2- JV / 2 ^| CT ), (4) 

3=1 <y 

where a J x is a Pauli matrix for j-th qubit. Consider the instantaneous eigenstates of H(t) with 
eigenvalues Afc(r) arranged in nondecreasing order at any value of r G (0, 1) 

H(r)\Mr)) = h(r)\Mr)), (5) 

here k = 0, 1, 2, . . . , 2 N — 1. Provided the value of T (the runtime of the algorithm) is large 
enough and there is a finite gap for all r G (0, 1) between the ground and excited state energies, 
Ai(r) — Ao(t) > 0, the quantum evolution is adiabatic and the state of the system |^(r)) stays 
close to an instantaneous ground state, |0o( r )) (up to a phase factor). The state 1 O ( 1 ) ) coincides 
with the ground state of the problem Hamiltonian Hp and, therefore, a measurement performed 
on the quantum computer at the final moment t = T (r = 1) will yield one of the solutions of 
COP with large probability. 

The standard criterion for adiabatic evolution is usually formulated in terms of minimum exci- 
tation gap between the ground and first exited states 

T » "XT2 — ' AA min = max [Ai(r) - A (t)] . (6) 

AA^ in 0<r<l 

Here the quantity £ is less than the largest eigenvalue of the operator Hp — Hp Jl^l and scales 
polynomially with iV in the problems we consider. 
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III. QUASICLASSICAL APPROXIMATION AND COMBINATORIAL LANDSCAPES 

In the computational basis © we have 

h = tY,eM (*\ -(i-t)J2 s M*, °-')> i] W) (o-'l> (7) 

<T <T,<T' 

here S[m,n] denotes the Kronecker delta-symbol and the summation is over the pairs of spin 
configurations cr and cr' that differ by the orientation of a single spin, d(cr, cr')=l, where 

1 - 

d(<T,<T') = -Y J Wj-°' J l (8) 

i=i 

denotes a so-called Hamming distance between the spin configurations cr and cr', that is the num- 
ber of spins with opposite orientations. Eq. © in the computational basis takes form 

A(r)0 CT (r) = tE v 4> v (t) - (1 - 7)^*^/), 1] 0^(r) (9) 

cr' 

(here we drop the subscript indicating the number of a quantum state in A and </>«,.) . In what follows 
we assume that typical energies E„ = O(N), but the change in the energy after a single spin flip is 
0(1). This assumption about the energy landscape holds for instances of the Satisfiability problem 
with the clause-to-variable ratio M/N = 0(1), the case of most interest for us (see the discussion 
in Sec. HI). 

We now consider a set of functions {Xi = C/(<r,Z), I = I,... ,/C}, referred to as (com- 
binatorial) landscapes, that depend on a problem instance X and project a spin configuration cr 
onto a vector {X{\ with integer- valued components. Prior to considering a specific COP here we 
make certain assumptions about the properties of landscapes and apply them to the analysis of the 
minimum gap in the QAA. 

In particular, we assume that, similar to energy, landscapes \X\ = Ci(cr,I)} are macroscopic 
functions, so that the typical values of Xi are 0(N), and possess a certain universality property 
in the asymptotic limit N — » 00. Specifically, the joint distribution of {Q(<x,Z)} over the spin 
configurations cr forming the 1 -spin-flip neighborhood of an "ancestor" configuration cr' depends 
on a problem instance X and spin configuration cr' only via the set of parameters {X[ = Cj (cr', X) } . 
We then define a quantity 

P({X l }\{X' l }) = ^ Yl rPt^a(<x,X)], (10) 

d(er,er')=l k=l 

x; = ch<x',x), 
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In effect, the above universality property of landscapes implies that the set of all possible 
spin configurations a is divided into "boxes" with coordinates {X/} where X/ = Ci(er), and 
P ({Xi} | {X/}) (TTOb represents the transition probability from box {X/} to box {X/}. In particu- 
lar, it obeys B ayes' rule 

p{{Xi] I W}) = p{{x[] | {x}) n({x,», (ii) 

where tt({Xi}) is the number of different spin configurations in the box {X{\. 
We consider energy to be a smooth function of landscapes 

E tr = E({X l }), X t = 0^,1), (12) 

so that \dE/dXi\ = 0(1). Furthermore, we assume that, on one hand, the change in Ci(cr,X) 
after flipping one spin is 0(1), for typical problem instances. On the other hand, we assume 
that correlation properties in a neighborhood of a box {Xi} described by P ({X{{ \ {X/}) vary 
smoothly with box coordinates on a scale 1 < \5Xi\ <C N. Therefore if we write the transition 
probability in the form 

P ({X/} | {X}) = p ({XI - X}; { Xl }) , { Xl = X/X}, (13) 

then p({ki}; {x;}) is a steep function of its first argument: it decays rapidly in the range 1 < 
\ki\ <C X for each /-component. However this is a smooth function of its second argument: it 
varies slightly when coordinates xi change on a scale \5x\\ 1. 

One can show that under the above assumptions the quantum amplitudes CT corresponding to 
the smallest eigenvalue depend on the spin configuration cr only via the coordinates of this box 
{Xi} to which it belongs. Then we look for the solution of © in the following form: 

= v( {Xi},r) Ci(er,l)}. (14) 

where |<£>({X}, r)| 2 gives the probability of finding the system in the box {X}. Plugging (TBI) 
into © and making use of (fTTl) .(fT2ll we obtain: 

A(rV(X, r) = tE(X)v(X, t) - (1 - r)N ]T L(X, X'}) V (X, r), (15) 

X' 

X={X 1 ,X 2 ,...,X k; }, (16) 
(hereafter we use the above shorthand notation for the set of landscapes). In (fT5t we introduced 



L(X,X) = L(X,X) = P(X\X)J^-, (17) 



PX 



p(x) = 2- N n(x) 
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where P(X) is a probability that a randomly sampled configuration er belongs to a box X. We 
shall look for a solution of (1131) in the WKB-like form 

^(X,r) = exp(-^(X,r)), (18) 

so that 

A(r) = r£(X) - (1 - r)iV^L(X,X> w(x ' r) - w(x '' r) . (19) 

X' 

We now introduce scaled variables (cf. (fT3t ) 



and also 



w(x,r) = i^(X,r), e(x) = if?(X), s(x) = 1 logfi(X), (21) 



where s(x) is an entropy function. Based on (fTTt and the properties of the transition probability 
(see Eq. (fT3l) and discussion after it) we assume that the sum over X' in ( fT9l is dominated by terms 
with |X' — X| = 0{1). Then we can use an approximation 

W(X', t) - W(X, r)^Vw (X - X) + 0(1/N), (22) 

where Vw = chv(x, r)/<9x. Plugging d22l) into ( fl9l and making use of Eqs. (fT3l.([T7l).(l2"Ot and 
d2TT) we obtain after some transformations: 

</ = /i(x,W;r), (23) 
fc(x,p;r) = e(x)-r^p(k;x) e - k ^/ 2+ P). 

k 

(here V* = <9s(x)/<9x). This is a Hamilton-Jacobi equation for an auxiliary mechanical system 
with coordinates x, momenta p = Vw, action w, Hamiltonian function h(x, p; T) and energy g. 
Using the symmetry relation 

p(k; x)e- k ' Vs / 2 = p(-k; x)e k Vs/2 , (24) 

that follows directly from Eqs. (fTTT) and (fTTt we obtain that the minimum of w(x, T) over x where 
Vw = necessary corresponds to the minimum of the functional: 



/( x ,r) = £ (x)-rY(x), 



(25) 
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where /(x, T) = h(x, 0, T) and 

£(x)=p(V S /2;x), p(y;x) = ^p(k;x)e- k ^. (26) 

k 

The summation in d23t and d26t is over components A;; of k in the range ki E (—00,00). In what 
follows, we shall refer to p(y; x) in (l26l) as a "Laplace transform" of p(k; x). 

We note that £(x) = L(X', X) and one can use Bayes rule and inequality of Cauchy- 
Bunyakovsky in (ITvb to show that that the positive-valued function i(x) is bounded from above, 
< £ (x) < 1. This shows that the analysis of the effective potential based on the WKB approxi- 
mation d22l is self-consistent in the asymptotic limit N — > 00. 

It follows from the above analysis that the ground-state wavefunction ^(x, T) = y?(X, r) 
is concentrated in x-space near the bottom of the "effective potential" given by the functional 
/(x, T), i.e. near the point x* (r) where /(x, T) reaches its minimum. In this region S ~ |x T A x, 
where matrix A is positive definite, and according to (TTHl) . the wavefunction has a Gaussian form 
with the width oc 1 /y/N. 

The ground-state energy g = g(Y) is given by the value of the effective potential / (l25t at its 
minimum 

<7(r) = /(x*(r),r), (27) 
9/(x,r)/ax| =o, /(x,r)>^(r). 

We note that as T — >• the shape of the effective potential /(x, T) approaches that of the energy 
function e(x) and therefore its minimum x*(r) — > x where x is a minimum of e(x). It can be 
shown that in this limit the ground-state eigenvalue approaches the minimum energy value e(x ) 
and the eigenvalues of A -1 approach zero (and so does the characteristic width of the wavepacket 
V>(x, T)). The spin configurations that belong to a box x in x-space correspond to the solutions of 
the optimization problem at hand. It is clear that one of the solutions can be recovered with high 
probability after a measurement is performed at the end of the "quantum annealing" procedure. 

Variational Ansatz: For cases in which the set of macroscopic variables {X[} is not suf- 
ficient (in statistical sense (fT5l ) to describe the dynamics of the quantum algorithm, one can 
still implement the above procedure as an approximation, using a variational method. Intro- 
ducing a Lagrangian multiplier A, one looks for the minimum of the functional F(cp, A) = 
(4>\H\4>) — A((0|0) — 1), using a variational ansatz (TBI for the wavefunction. The solution of 
the variational problem is provided by Eqs. dT8b-d27T). The smallest eigenvalue g (l27t corresponds 
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to the value of the Lagrange multiplier at the extremum, A = rNg, and the maximum of the 
variational wavef unction corresponds to the minimum of the effective potential / d25l) . 

A. Global bifurcations of the effective potential 

However, in the case of a global bifurcation where the effective potential f(x, V) possesses 
degenerate or nearly degenerate global minima, the answer is modified. If for some value of 
T = r*, a global bifurcation occurs, in our example this would mean that for this value of T, two 
values of x, x+ and x~ give a global minimum to /(x, T). In such a case, the smallest eigenvalue 
is not doubly degenerate; rather an exponentially small gap AA min between the ground and first 
excited state is developed, itself being proportional to the overlap between two wave-functions, 
peaked around x+ and x~ respectively. 

To estimate the overlap we note that at the two global minima of the effective potential 
/(x, r*) correspond to the two coexisting fixed points of the Hamiltonian function in d23l with 
zero momentum and the same values of energy g, 

df/dx = dh/dx = dh/dp = (28) 
x = x^, p = pf = 0, g(x,p;T*) = gt = g~. (29) 

Then to logarithmic accuracy we have 

1 f°° 

-logA to = / dt'[x(t')p(t')-h(x(t'),p(t'))} + 0(l/N), (30) 

J — oo 

where ( x(t), p(t) ) is a heteroclinic trajectory connecting the two fixed points of (l23l 

x(t) = dh/dp, p(t) = -dh/dx, (31) 
x(t -> ±oo) = xf, p(t -> ±oo) = 0. 

From the algorithmic perspective this means that when T gets close to I" 1 *, it has to change 
exponentially slowly (cf. Sec. m and Eq. ©). This could be called a critical slowing down in 
the vicinity of a quantum phase transition. If simulated annealing (SA) is used and a similar 
phenomenon occurs, the value of the temperature T* is the point where a global bifurcation occurs 
in the free energy functional 

f(x,T) = e(x)-Ts(x). (32) 
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By comparing the free energy functional (l32l) with the functional (l25t corresponding to "quantum 
annealing" (QA), we note that in QA the quantities T and £(x) play the roles of temperature and 
entropy in (SA), respectively. 

We note in passing that a similar picture for the onset of global bifurcation that can lead to 

nn 

the failure of QA and (or) SA was proposed in [15, 16] for the case where the energy E a is a 
non-monotonic function of a single landscape parameter, a total spin Ylj=i &j- 



T n this case the 



dynamics of QA can be described in terms of one-dimensional effective potential y/7 



18] 



IV. THE MODELS 

An instance of a Satisfiability problem with N binary variables committed to M = •yN con- 
straints (where each constraint is a clause involving K variables) can be defined by the specifi- 
cation of the following two objects. One of them is an M x iV matrix Q, the rows of the matrix 
are independent /^-tuples of distinct bit indexes sampled from the interval (1, N). The m-th row 
of Q defines the subset of the K binary variables involved in the m-th clause. The second object 
is a set of boolean functions B = {b m }, with each function encoding a corresponding constraint. 
A function b m = b m [a g i ,a g 2 , . . . , a g ] is defined over the set of 2 K possible assignments 
of the string of K binary variables involved in the m-th clause. The function returns value 1 for 
assignments of binary variables that satisfy the constraint and for bit assignments that violate it. 
Then the energy function equals to the number of violated constraints 

M 

E (T = E (T {l)=M-y^b m [a g a Q 2 ,...,a g ], (33) 

m=0 

here X = (Q, B) denotes an instance of a problem. 

The matrix Q defines a hypergraph Q that is made up of the set of iV vertices (corresponding 
to the variables in the problem) and a set of M hyperedges (corresponding to the constraints of 
the problem), each one connecting K vertices. An ensemble of disorder configurations of the 
hypergraph corresponds to all the possible ways one can place M = "fN hyperedges among iV 
vertices where each hyperedges carries K vertices. Under the uniformity ansatz all configurations 
of disorder are sampled with equal probabilities (i.e., rows of the matrix Q are independently and 
uniformly sampled in the (l,N) interval). 



Boolean functions b m may also 



being random K-SAT problem 



OH 



3e generated at random for each constraint with an example 



However here we consider slightly different versions of 
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the random Satisfiability problem that are still defined on a random hypergraph Q but have a non- 
random boolean function b m = b, identical for all the clauses in a problem. One of the problems 
is Positive 1-in-K Sat in which a constrain is satisfied if and only if exactly one bit is equal 1 and 
the other K-l bits are equal 0. The boolean function b for this problem takes the form 

K 



b[«i, «2, • • • , olk] = 5 
±1, p=l,2,. 



E 1 



a 



P =i 
.,K. 



(Positive 1-in-K SAT). 



(34) 



We shall also consider another problem, Positive K-NAE-Sat, in which a clause is satisfied unless 
all variables that appear in a clause are equal ("K- Not- All-Equal-Sat"). The boolean function b 
for this problem takes the form 

" K 

b[«i,a 2 , . . . ,a K ] = 1 - ^ 5 



8=±1 



j2 1 + sap ,o 



, P =i 



(Positive K-NAE-S AT). (35) 



Both problems are NP-complete (Appendix©. It will be shown below that they are characterized 
by the same set of landscape functions. 



V. LANDSCAPES : ANNEALING APPROXIMATION 



For a particular spin (cr) and disorder (£?)) configurations, all clauses can be divided into 2 K 
distinct groups according to the values of the binary variables that appear in a clause. We will label 
the different types of clauses by vectorial index a = {a±, . . . , ax}, ot p = ±1. We now divide 
the set of 2^ spin configurations into boxes identified by certain numbers of clauses of each type, 
N A4 a , and also by the Ising spin in a configuration Nq 

M K 



M a = M a ((r,G) 



ton 



a P 



m=l p=l 



1 N 



(36) 



(37) 



3=1 



Different boxes correspond to macroscopic states defined by the set of parameters (q, {Ai a }) with 
q E (—1, 1) and ^2 a M- a — 7- The energy function can be expressed via d3*6l) as follows (cf. 



K 



K 



a, 



(38) 
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where the form of the coefficients ( m depends on the problem: 

6[m, 1] (Positive 1-in-K SAT) 



Cn 



(39) 

5[m, 0] - 5[m, K] (Positive K-NAE-SAT) 



In the following we compute an approximation to the effective potential (l25t . using the land- 
scape functions (l36b . (l37t . According to (1261 it depends on the entropy function s(q, {M. a }) and 
the transition probability dT3l) between different macroscopic states. Recalling that variables q and 
Ai a are normalized by the factor N we study the probability of transition, p(n, {r a }; q, {Ai a }), 
from the state (g, {M a }) to the state (q + n/N, {M a + r a /N}). The Laplace transform of p 
with respect to n, {r a } has the form (cf. (l26l) ) 

p(0,{2/a};g,{^«})= e-^^^jrJi^MJ), (40) 

n, {r a } 

We assume that all binary variables are also subdivided into distinct groups based on their value 
a = ±1 and a vector k with integer coefficients k v a indicating the number of times a variable 
appears in a clause of type a. in position p. Clearly, consistency requires that k p a = unless 
a p = a. We now define a quantity c CT) k which is equal to the fraction of spins with given a, k. For a 
spin configuration cr there exists a set of coefficients {c a ,k} with elements of the set corresponding 
to all possible values of a and k (there will be many 0's in a set for each spin configuration). 
In general, there are exponentially many sets {c a ^} that correspond to a macroscopic state (q, 
{M a }) 

^2ac aM = q, ^2k p a c aM = M a (p = 0,l,...,K). (41) 

cr,k cr,k 

Coefficients {c CTj k} are concentrations of spin variables with different types of "neighborhoods". 
We shall assume that in the limit of large the distribution of coefficients c CT k corresponding to 
the same macroscopic state d4TT) is sharply peaked around their mean values (with the width of the 
distribution oc N~ 1/2 ). 

Under the above assumption we can immediately compute the Laplace-transformed transition 
probability (l40l in terms of the coefficients c CT k- Indeed, consider flipping a spin with value a 
and neighborhood type given by vector k. This will change the total spin by —2a and for each 
clause of type a. and index p E (1, K) the value of NM. a will decrease by . On the other 
hand, for the clause type a' = a(p, a.) obtained by flipping a bit in p-th position in ck, NAi a > is 
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correspondingly increased by k v a . Hence the Laplace-transformed transition probability is 

P (0, {y<*}; q, {Ma}) = ^2 c<T > k exp 

cr,k 

where the coefficients c CT) k are set to their mean values in a macroscopic state 

A. Entropy and coefficients c CT; k in a macroscopic state defined by q and {M a } 

Here we use the annealing approximation to estimate the mean values of c^k and also of a 
macroscopic state (q, Ai a )- We start by introducing the concept of annealed entropy. Let J\f be 
the number of spin configurations subject to some constraints. In general, it is a function of the 
disorder realization. The annealed entropy is defined as the logarithm of its disorder average: 
s ann = In (A/). Note that for the correct, quenched, entropy the order of taking a logarithm and 
disorder average is reversed. 

Since in the random hypergraph model all disorder configurations are equally probable, an- 
nealed entropy is given as s ann = lnJV^g — InA/g, where Ms,g is the total number of spin and 
disorder configurations and J\fg is the number of disorder configurations. 

For enumerating all possible disorder configurations we depart slightly from the traditional 
random hypergraph model. In our model all clauses are ordered (two disorder configurations where 
any two clauses are permuted are deemed different), clauses can be repeated (the same clause can 
appear twice), the order of variables in a clause is important (two disorder configurations are 
different if the order of variables in any clause is changed), and finally, variables can be repeated 
in a single clause. This change does not alter the underlying physics, since the probability that 
two identical clauses appear is infinitesimal, and a variable enters a clause twice in at most 0(1) 
clauses, which can be safely neglected. As regards the distinction between the disorders with 
permuted clauses, this only introduces a combinatorial factor which cancels out. The advantage is 
that each disorder can be represented as a sequence of M i^-tuples of integers from 1 to N. 

We will first compute the annealed entropy of a macroscopic state (q, {Ai a }) under addi- 
tional constraints: we fix the values c CT; k and compute the annealed entropy as a function of 
q, {-Ma}, {ccr,k}- Recalling that Ai a are the numbers of clauses of a given type scaled by N, 
and the total number of clauses is jN, we obtain the number of joint spin-disorder configurations 
as a product of the following factors: 

(i) the number of ways to assign types to clauses (N^)\/ Y[ a (NM a )\, 



a yot(p,ot) 



(42) 



p.ot 
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(ii) the number of ways to assign types to variables AH / k(-^ c o-,k) !, 

(iii) for all p, a, the number of ways to permute the appearance of variables in p-th position of 
clauses of type a: {NM a )\/ Y[ a , v {K}) Nc ^, 

Consequently, the annealed entropy is given by 



Sarin [{c CTjk }; Q, {M a }} = - ^ C CT)k In 



cr,k 



C CT ,k J) 



+ l)J^>l ce ln>l ct + 7ln7-7ir. 

(43) 

In the large N limit we replace c a ^_ by their annealed averages, i.e., the values that maximize the 
annealed entropy. In its simplest form, we place no constraints on c CTj k except consistency require- 
ments d4TT) . Associating Lagrange multipliers A and ln/x^ with these constraints, the expression 
for the entropy can be rewritten as 

Sannfe, {M a }] = min \ -Xq + y2M a \n^- + \nZ[\, {fj, p a }) 

A, a* { p a Vex 

- M a \nM a + 7ln7 - 7K (44) 
The values of c CTi k are given by 



Z 

p, OL 



and Z is given by 



\ at p / \ a p / 



The values of the Lagrange multipliers A, \i v a are related to q, {M a } via 

dlnZ 



d\nZ 



q, (47) 



From here we obtain the expression for the Lagrange multiplier /x£ 

M a _ 1 + a p q 
1K~ 2 ' 

Then introducing a new notation 



El±a p 
— — /C 



2 

p,a 



(49) 
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we obtain 



Z = e V+ + e~V- , a± = (51) 

1 ± q 



Then the entropy can be rewritten in the following form 



s a nn[q,{M a }] = -Xq+M+te^+M-te^+lnZ-^MnkiMn+^-ryK. (52) 



We now use the following equations 

_-A_U_ ry^-q 



e V+ = Z—^, e'V" = Z ^ (53) 

2 2 

and obtain the expression for the second Lagrange multiplier A 

—Xq = -±±* In l -±l - ?-± In - In Z + 7 tf. (54) 

y 2 2 2 2 

Upon substitution of A from the above into the expression for s ann (1521 we finally obtain the 
annealed entropy 



Sann[g, {^a}] = - 9 tami" 1 q - In ^ q + A4+ In + A<_ In 

- M a \nM a + 7 In 7. (55) 

Also the coefficients c CTik are given by (I431).d4*6l) with Lagrange multipliers given in (14*91) and (1541. 



B. Effective potential 

Consider a factor £(x) = (9/ / <9r) e d25t . (EST) in the expression d25l) for effective potential with 
x = (q, {AA a }). It follows from (|2"61) that to find this factor we need to evaluate the Laplace- 
transformed probability (1401421) ) at 

= ^ds ann /dq, y a = ^ds ann /dM a . (56) 

This is where the Lagrange multipliers come in handy as we can immediately claim that 

ds ami /dq = -A (57) 
d Sam /dM a = Vln-^-lnM,. (58) 
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Note that in differentiating with respect to M. a above we omitted the constant term. This is 
permissible since only differences dq ann /dM a — dq ann /dM a ' appear in Eq. (l42b . A further 
refinement is to write 



E ln ^x = E ln l - ± r 1 = K ln + 5> tanh " (59) 

p ^ a P V 
Using this in the Eqs. (I26l) . (l42l . we obtain 



Since | J2 p >( a P' ~ a ' P ') = °p (where ex' is obtained from ex by flipping p-th bit) and also 



MJ& = y/1 2 t- e"'^-**, (61) 
the expression is considerably simplified 



{M a }) = |exp (— == VM a M a > J 



where the sum is over pairs (a, ex') that differ in exactly one position 

K 

2 



P =i 

To evaluate Z we write 



2 f M+ M. 



(62) 



1 K 

-^|a,-a;.| = l. (63) 



Z = Ve^+e^- = ; exp — H (64) 

yr 3 ^ a/w^ + 9 1-?/ 

and the expression for i becomes 

*(?, {M a }) = exp ^ j== _ - —j . (65) 

here M.± are given in (l50l) . 

We note that the effective potential f(q, {M a }) = e({M. a })— T£(q, {Ai a }) is symmetric with 
respect to permutation of individual components in {Ai a } corresponding to different orders of - 
l's and +l's in the vectorial index ex. We look for the minimum of f(q, {M. a }) using symmetric 
ansatz 

M a =( K ) M m , m = ^i_^. (66) 

\ m J p=1 
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where m is the number of -l's in ex. Substituting (l66t into (l65l) and rewriting 




K 1 + qJ2* =0 (K-2m)M^ 



(67) 



1 -g 2 

where we defined £(q, {M rn }) = {.M a }). The effective potential is then 

f(q, {M m }) = e ({M m }) - T£(q, {M m }) (QA), (68) 

with energy given in d38t . In the case of the SA algorithm the corresponding free-energy functional 
(Ellis 

f(q, {M m }) = e({M m }) - Ts(q, {M m }) (SA), (69) 
where the entropy function equals 



/ 1 

s(q,{M m }) = - gtanh- 1 g+( 7 K-l)ln V ~ g 

X \ K M 

J2( K ~ 1rn)M m j tanh" 1 q - ^ M m In (70) 

\m=0 / m=0 \m/ 

If we were to use an even smaller set of macroscopic parameters (e.g. only the energy e) 
we can still employ formula d67t with the proviso that unspecified variables should be taken to 
equal their most likely values, i.e. those that maximize the entropy s(q, {M m }) not the landscape 
{M. m }). For example, in the case of energy-only landscapes, I = 1(e), the values q, {M m } 



that maximize s(q, {M m }) for a given energy e and number of hyperedges 7iV (J2m=o M m = 7) 
should be computed and then substituted into the expression for i doTt . 

We compute, within the annealing approximation, the point of static transition 7 C (cf. Fig HI), 
where the entropy of the macroscopic state with zero energy vanishes, s(0) = 0, and the dynamic 
transition 7^; for connectivities 7 > 7^ an effective potential (l68l) exhibits a global bifurcation for 
some r = r*. The resulting values are given in Table0](see also Figs. E]and[21). Note that in l-in-3 
SAT and K-NAE-SAT for (K=3,4,5) we find no dynamical phase transition before the satisfiability 
threshold (cf. Fig. [T]). 

In Fig.|4]we plot time variations of the landscape parameters, M rn = M* m , corresponding to the 
global minimum of the effective potential. In Fig. |5]we plot a time-variation of the scaled ground- 
state energy g given by the value of the effective potential at its minimum. Singular behavior 
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K 


3 


4 


5 


6 


7 


8 


9 


10 


1-in-K 


Id 
7c 


0.805 


0.650 
0.676 


0.557 
0.609 


0.475 
0.548 


0.416 
0.500 


0.371 
0.461 


0.335 
0.428 


0.305 
0.400 


K-NAE 


Id 
7c 


2.41 


5.19 


10.7 


19.8 
21.8 


34.9 
44.0 


61.7 
88.4 


109 
177 


196 
355 



TABLE I: Annealing bounds for dynamic (7^) and static {y c ) transition for positive 1-in-K SAT and positive 
K-NAE SAT for different values of the number of variables in a clause K. No value (-) indicates the absence 
of dynamical transition. 
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FIG. 2: Static j c (circles) and dynamic 7^ (di- 
amonds) transition for positive 1-in-K SAT for 
various values of K. 



FIG. 3: Static j c (circles) and dynamic 7^ (di- 
amonds) transition for positive K-NAE SAT for 
various values of K. 



corresponding to the first-order quantum phase transition at certain r = t* (r = T*) can be clearly 
seen from the figures. Plots in Figs.|4|and|5]correspond to precisely the static transition 7 = j c for 
the case of K = 4 in 1-in-K SAT problem. In the region 7d < 7 < 7 C there are exponential (in N) 
number of solutions to Satisfiability problem but the runtime of the quantum adiabatic algorithm 
to find any of them also scales exponentially with N. This is a hard region for this algorithm. We 
note, that in the limit of K — > 00 the annealing approximation becomes exact. Together with the 
fact that for large K ^ d and 7 C seem to be distinctly different provides evidence that this result 
(existence of hard region for quantum adiabatic algorithm) is robust. 
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FIG. 4: Plots of the landscape parameters M m = M* m at the global minimum of the effective potential, vs 
t for K = 4 and 7 = 7 C (1-in-K SAT problem). Curves labeled 0-4 correspond to M*o/7 through M^/^y. 
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FIG. 5: Scaled energy of adiabatic ground state g$ vs r 
for K=4 and 7 = 7 C (1-in-K SAT problem). 

VI. UNIVERSALITY PROPERTY FOR TRANSITION PROBABILITIES 

Here we study the universal features of the transition probability in (flOb for the set of macro- 
scopic variables corresponding to the (normalized) total Ising spin q and numbers of clauses of 
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FIG. 6: Results of numerical simulations and their comparison with theory. Depicted are Laplace trans- 
forms of M.i for l-in-3 SAT. Numerical results: curves that have different colors correspond to different 
random problem instances; curves of same color correspond to different random bit strings. The dashed 
black line is a theoretical result based on the annealing approximation. The insets (a)-(e) depict instances 
with 10 3 , 10 4 , 10 5 , 10 6 , and 10 7 binary variables. Since the error is not recognizable we replot in (f) a 
magnified section of inset (e). The bit strings were sampled with q = 0.422, Mq = 0.048, Mi = 0.416, 
Mi = 0.123, M3 = 0.013, corresponding to M/N = 0.6. These values correspond to the energy 
and they are shifted by 10% from the most likely values of q, {M m } for this energy (this shift is > N 1 / 2 ). 
We also note that for l-in-3 SAT numerical simulations give static phase transition at 7 C 0.63). 
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different types {M rn } d38t (the type of a clause is equal to the number of unit bits involved in the 
clause). For simplicity, we shall focus in this section on the case K = 3 only. 

To clarify the above choice of macroscopic variables we consider an auxiliary quantity: a 
conditional probability distribution of the macroscopic variables (q, {M m }) over the set of all 
possible configurations cr obtained by flipping r bits of the configuration cr'. The first moments of 
this distribution corresponding to M m , 



(^) J^Ko-'-cr), r] M m {(T,l), m = 0,...,K, (71) 



can be easily computed by counting the number of ways one can flip r bits in configuration cr' to 
transform a i^-bit clause of m! type (i.e., with m! unit bits) into a clause of the m-th type 

(N\~ l A _ fm'\f K-rri \( N-K \ 
/i m = V M m , , (72) 

V r / i q \P / \ m ~ m + P/ \r — 2p — m + m J 

(here we use the convention ( n ) = for m < and m > n). In the double sum above values of 
Aim' are multiplied by the number of possible ways to flip three groups of bits: p unit bits in a 
clause of m'-type, p + m — m' zero bits of this clause, and r — 2p — m + m! bits of the configuration 
cr' that do not belong to the clause. Similarly, one can show that the first moment corresponding to 
the variable q equals q'(l — 2r/N). It is clear that dependence of the first moments on the ancestor 
configuration cr' is only via the variables q', M' m for that configuration. 

In the limit, r > 1, the above conditional distribution has a Gaussian form with respect to q 
and M m . Elements of the co variance matrix E^' 9 '(cr') = 0(r), and correspondingly, the charac- 
teristic width of the distribution is 0{r x l 2 ). For a configuration cr' randomly sampled in the box 
(q, {M m }) the r.m.s. deviation of the elements of S^;' 9 '(cr') from their mean values in the box is 
0(N 1 ^ 2 ). It is clear that in the limit r ^> N 1 / 2 the covariance matrix elements can be replaced 
by their mean values for the macroscopic state (q, {M m }). Therefore in this limit the conditional 
distribution after r spin flips starting from some macroscopic state depends only on the values of 
(q, {M m }) in this state (universality property). 

One can show that for r N 1 / 2 the conditional distribution after r spin flips can be expressed 
via the distribution (fTDI) with r = 1, using a standard convolution rule. Implicit in our derivation 
of landscapes is the relation between r = 1 landscapes and a set of quantities {c CT k}. Universality 
of landscapes should be interpreted as the fact that {c CTj k} are self-averaging. Had we included 
only the energy instead of a full set of parameters, we would not have expected to see such self- 
averaging in the so-called replica- symmetry -broken phase. It is possible that inclusion of the full 
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set of landscape parameters assures universality. We performed a series of numerical studies to test 
this hypothesis. In Figure |6] we present the results of numerical simulations and the comparison 
with analytic results within the annealing approximation. One can see that the property of self- 
averaging holds and that even annealing approximation provides very good description. 

VII. RANDOM GRAPH ENSEMBLES WITH REDUCED FLUCTUATIONS 

One can easily point out to a major deficiency of the annealing approximation - the fact that it 
fails to correctly predict the satisfiability transition. While part of it can be attributed to the fact that 
real entropy is slightly different from what is predicted by annealing approximation, major source 
of error is the incorrect assumption that the entropy vanishes at the satisfiability transition. That 
the entropy does not vanish can be easily seen by examining the structure of the random graph. At 
any finite connectivity (above percolation) it consists of one giant component and a large 0(N) 
number of small 0(1) components. Each small component contributes 0(1) contribution to the 
entropy for a total of 0(N) contribution, hence the entropy is in fact positive at all connectivities, 
including the satisfiability threshold. 

A. Concept of a core 

An improvement over annealing bound for satisfiability threshold is possible. Note that clauses 
outside of giant component do not affect satisfiability and hence can be disregarded. Similarly, one 
can identify irrelevant clauses and remove them. Irrelevant clauses are those that can be eliminated 
without changing the satisfiability of the entire problem. We shall illustrate the identification of 
irrelevant clauses based on local properties for positive i-T-NAE SAT and positive l-'m-K SAT. 

1. Positive K-NAE SAT 

For i^-NAE SAT we can identify variables that do not enter any clause and remove them 
without affecting the satisfiability of the formula. Any variable that appear in exactly one clause 
can also be eliminated together with such clause, since the value of that variable can be adjusted 
to satisfy that clause. As one removes such clause, other variables may become candidates for 
removal. One can write an algorithm that iteratively removes variables that appear in no or exactly 
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one clauses until no such variables remain (see Fig. 13). In fact using such algorithm improves 
running time of classical algorithm. The cost of this additional preprocessing is negligible: by 
using special data structures this algorithm can be made to run in 0(N\n N) time. 

It is not surprising that the result of running such trimming algorithm is a set of clauses and 
variables with a condition that every variable appears in at least two clauses. However, the statisti- 
cal properties of the remaining core are not immediately apparent. As a first step, observe that the 
resulting core is independent of the order in which clauses and variables were removed. Hence we 
can study one specific algorithm. 




FIG. 7: Example of trimming algorithm for 3-NAE SAT. Variables are represented graphically as vertices 
and clauses are represented as triangles. Incomplete triangles represent connections to the remainder of the 
graph (not shown). Shaded clauses are removed by the trimming algorithm. 

First, we eliminate all variables of degree (that appear in no clauses). In the initial problem 
all instances with A" variables and M clauses were equally probable. The number of vertices 
of degree 0, A" does fluctuate, but for every specific value of N all instances with N — N 
variables, and M clauses, and the additional property that each variable appear in at least one 
clause, are equiprobable. This property is referred to as uniform randomness. As a next step, 
identify variables of degree 1 (appearing in exactly one clause). Although their number Ai also 
fluctuates, for any fixed N\, N and M all instances with A" variables of which N\ variables have 
degree 1 and M clauses are equiprobable. 

At each subsequent step of the algorithm we choose at random a variable of degree 1 among 
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Ni candidates, and delete it together with the clause in which it appears. If the degree of any other 
variable in that clause becomes 0, it is also deleted. Some variables which had degree 2 previously 
may become degree 1 variables. Obviously M' = M — 1, but the values of N' and N[ cannot 
be predicted. But since the variable chosen for deletion was chosen at random (among all degree 
1 variables), one can show that although the values of N' and N[ cannot be predicted, for every 
fixed set of {N', N[, M'} all instances are equiprobable. 

At the end of the algorithm Ni = 0, hence for fixed N', M', all instances with N' variables, M' 
clauses and the condition that each variable appear in at least two clauses are equiprobable. Since 
average changes in N, Ni and M at each step were 0(1) and corresponding deviations were also 
0(1), after O(N) steps needed for completion of the algorithm, fluctuations in resulting N' and 
M' are only 0(yf~N) by central limit theorem and can be neglected in our analysis. 

There are several methods to derive the average resulting values of N' , M' . The most straight- 
forward is to study the evolution of N, N 1 and M by solving differential equations for their average 
values. This is somewhat tedious and less general. We instead resort to a different method of self- 
consistency equations. For every instance one can identify a set C of variables that form a core, 
i.e. those that will not be delete by the trimming algorithm. We now extend that set to C using the 
following rules: 

1 . If a variable belongs to C, it also belong to C. 

2. If K — 1 variables in some clause belong to C, then the remaining variable must also belong 
toC. 

The minimal set C is obtained from C by iteratively adding variables according to these rules. Let 
the number of variables in set C be equal to qN with q < 1. Fluctuations in q are 0(1/ y/N) and, 
consequently, are neglected. 

Now introduce (N + l)-st variable and compare systems with N variables and N+l variables. 
Together with (N + l)-st variable introduce AM clauses each involving that variable and (K — 1) 
random variables. The number AM is Poisson random variable with parameter Kj. We can 
compute the probability that the new variable belongs to C of the enlarged instance. The new 
variable belongs to C if for at least one of AM clauses, K — 1 variables other than the new 
variables are not in C at the same time. Since these clauses are connected to random vertices, 
that probability can be expressed via q alone and, owing to the fact that AM is Poisson, equals 
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1 — exp(—K r yq K 1 ). Self-consistency requires that this probability be equal to q: 

l-q = e - K ^ K '\ (73) 

For practical purposes we must seek the largest solution to this equation. 

Note that the new variable belongs to the core (set C) if not for one, but for at least two of AM 
clauses, K — 1 variables other than the new one are in C. The probability of that 

p = 1 - e-^- 1 - K iq K - x e- K ^ K ~ X =q- K iq K -\l - q) (74) 

is equal to the average value of N'/N, and the average connectivity of the new vertex (under the 
condition that it is in C) enables to compute M'/N': 



1Q K 



M' IN' = - — = — ±— (75) 

(better expressed as M'/N = ^K^q K ~ x (l - e - K ^ qK ^ = -fq K ). 

We must mention that the core for i^-NAE SAT is in fact identical to that for K-XOR SAT 
The latter has been analyzed by a different method J^lll . 



2. Positive 1-in-K SAT 



The situation is slightly more interesting in case of positive l-'m-K SAT Setting a variable that 
appears in only one clause does not guarantee that the clause will be satisfied, hence the clause 
cannot be eliminated. Although it cannot be eliminated, it can still be simplified by eliminating the 
variable. The remaining instance will contain not only fT-clauses but also one (K — 1) -clause. The 
allowed combinations of variables for (K — 1) -clauses should be those that have sum of variables 
equal to either or 1. In that case the value of the deleted variable can always be adjusted so that 
the sum of variables in i^-clause is exactly 1. 

Similarly, at some point (K — l)-clauses can be converted to (K — 2)-clauses, (K — 3)-clauses, 
down to 2-clauses. In each case the allowed combinations shall have sum of spins either or 
1. Finally, if one can identify a variable that appears in 2-clauses exclusively, this variable can 
be eliminated together with corresponding 2-clauses, since setting it to shall satisfy all these 
clauses. As an illustration, see Fig. [8]which depicts K = 4. 

Once again, one can write trimming algorithm that eliminates such variables and clauses. As 
before, one can show that all remaining instances must have the property that each variable appear 
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FIG. 8: Example of trimming algorithm for l-in-4 SAT. Variables are represented graphically as vertices and 
4-, 3- and 2-clauses are represented as rhombi, triangles and edges correspondingly. Incomplete polygons 
represent connections to the remainder of the graph. The figure depicts evolution of part of the graph under 
trimming algorithm. 

in at least two clauses of any length and at least one clause of length greater than two; and that for 
fixed number of vertices N', number of 2-clauses M' 2 , (K — l)-clauses M' K _ X , fT-clauses M' K , 
all such instances are equiprobable. 



Without going into details (available in our upcoming paper H22I0 . we should mention that 
equations now involve two variables q and q', owing to separate treatment afforded to 2-clauses 
and clauses of length greater than 2. The self-consistency equations are: 

l-g' = e -K 1 [l-{l-q) K - 1 -{K-l)q{l-q) K -^\^ (76) 

l-q = (l_ g ') e -^(^-l)T 9 '(l-^- 2 ; (77) 



with implicit understanding that the largest solutions q, q' are sought. 
The expected values of N' as well as M' k are given by 

N'/N = q' - #7 [1 - (1 - qf- 1 - (K - l)q(l - q) K ~ 2 ] (1 - q), (78) 

M' 2 /N = ( K \ 1 q n {l-q) K -\ (79) 



M'JN = (^j iq k (l - q) K - k fork ^ 3. (80) 

Note that results for the core for positive K-NAE SAT and positive l-'m-K SAT allow for a com- 
pact formulation through introduction of generating function of allowed variable degrees. 

For the case of i<"-NAE SAT it is G{jj) = ^ d>2 l Jid I d\ = e M — 1 — \i, since degrees of all 
variables are at least 2. In terms of this G(n), equations can be written as 

N'/N = e^G(fi) (81) 

M'/N = —fie-"G'(fi) (82) 
K 
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with n = Kjq K 1 , and the equation on q is 

q = e-^G'i/j,). (83) 

For the case of 1-in-fT SAT, G is a function of (K — 1) variables reflecting the fact that the degree 
of each vertex is a vector (d 2 , ■ ■ ■ ., dx) where d k counts the number of appearances of said variable 
in /c-clauses. The generating function for the remaining core is 

K 

G(ji2, fi 3 ,...,fi K )= « - - IM- (84) 

i=3 

There appears a set {g 2 , . . . , q K } such that q k = (d/dfi k )G({^i}). Obviously, q 3 = q 4 = . . . = q K . 



Setting 



N' and M' k can be compactly written as 



& = 2( K \ iq \l-q) K -\ (85) 



tik = k\ ^\q k -\l ~ q) K ~ k folk > 3, (86) 



N'/N = e -2^ «(?({//«}), (87) 
M^/iV = I^e-Ef^w^^^.}). (88) 

B. Improved annealing bound 

Once irrelevant clauses have been removed we expect the entropy at the satisfiability transition 
to be much closer to zero. We have already made predictions for the number of remaining vari- 
ables and clauses after trimming algorithm starting from a random graph of connectivity 7. We can 
compute the entropy of the remaining core. Using the critical connectivity at which the annealed 
entropy of the core alone becomes as an estimate of satisfiability transition is an improvement 
over traditional annealing approximation. We are motivated by two contributing factors: (i) the 
rigorous mathematical proof that the disorder is relevant for satisfiability transition relies on pres- 
ence of irrelevant clauses, hence their removal can sometimes make disorder irrelevant, and (ii) 
applied to K-XOR SAT, this improved bound becomes exact 12 ill . 

Doing annealing approximation on a core is hampered by the fact that clauses can no longer be 
assumed uncorrected. However, the technique we introduced in this paper is well-suited for this 
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example. The annealed entropy equals 

*an„ = ^7 In Af(s, J) - ^ lnjV(J), (89) 

where jV( J) counts the number of possible disorders and J\f(s, J) is the total number of allowed 
spin-disorder configurations. We now separately consider K-NAE SAT and 1-in-K SAT. 



1. Positive K-NAE SAT 



Consider the number of possible disorder realizations with N' variables, M' clauses with the 
constraint that each variable appear in at least two clauses. Each variable can be characterized by 
the vector k with K components, each component k p counting the number of clauses where said 
variable appears in p-th position. Treating disorders that differ by permuting different clauses or 
permuting variables within a clause as different and fixing the fractions of vertices with particular 
realization of k to c k , the normalized logarithm of the count of possible disorders is 



-^lnAT(J) = S ;[iV',M';{c k }] = -$> k ln 



+ K^\nM'. (90) 



N 

We must maximize this expression with respect to c k subject to constraints 

J2h^ = M'/N' (91) 

k 

and that c k = for J2 P k p < 2. 

Introducing Lagrange multipliers /j, p , the expression for c k becomes 

Echoing the constraint that J2 P K ^ 2, the normalization factor is G(^2 where G(x) = 
e x — 1 — x - the familiar generating function of the core. Using dual transformation we can rewrite 
the entropy as 

s' [N', M'\ = min ^£(M'/N') In + }aG (g^ + ™L \ n N'. (93) 

Comparing this expression with Eq.() we observe that the minimum is achieved for \x v = j^n = 
r yq K ~ 1 , 7 being the connectivity of untrimmed random graph. 
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Now we proceed to count all possible combinations of spin assignments and disorder realiza- 
tions ftf(s, J). Corresponding entropy is in general a function of {.M^,} where M.' a is the number 
of clauses with variables having values described by a vector {a p } = a, divided by N' so that 
J2 a Ma. = M'/N'. In complete analogy with results of section IVAl we obtain 



s? 1 [N',M';{M , a };{c ir , k }] = -$> ff , k ln 



o-.k 



p,a 



K f^N: (94) 



N 

Optimizing this with respect to c CT; k subject to familiar constraints 

Ys kP c*c^ = M a (95) 

<r, k 

and the constraint that = if ^ a k v a < 2, the expression can be equivalently rewritten as 
s[[N',M';{M' a }] = mmlj2 M '^^ + ^ Z ^}}} + 

^ a (. p,ct ) ct 

+ ^N' (96) 
with Z[{[i p a }] given instead by 

\p,ct / \P,ct / 

Labeling arguments by fi + and /x_ respectively allows us to simplify this to 

( ML M' 

s ' {{M' a }) = min <^ M' In — ± + M'_ In — + ln[G(// + ) + G(fi. 
I /'• /' 

ct 

with^v = E a (E^K; i]) M' a and = £ Q (Epfc -i]) -K*- 

The correct annealed entropy is given by the difference of these expressions 

s' anQ [N', M'; {M' a }] = s[[N', M'\ {M' a }} - s' [N', M'\ (99) 

This enables us to find the improved bound on satisfiability threshold (where H ami = 0) or to 
compute the most likely values of {A4^} for particular energy E. 

Note that for positive fT-NAE SAT corrections to static threshold due to this improved approx- 
imation are minute, since the transition happens at large connectivities, where simple annealing 
approximation adequately describes the transition. 
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2. Positive 1-in-K SAT 

For l-in-.fr SAT the derivation is quite similar. An important addition is that the clauses can 
have any length from 2 to K, hence index k should reflect that fact. Since calculations are similar 
in spirit, we shall only provide the results 

s' [N>, {M' k }] = min L§ In *^ + In G( W) } + £ ^ In N>, (100) 
with G({fi k }) = exp (Ef= 2 /^) ~ e M2 - Ef= 3 J"fc> and 

K r i K kM> 

M' k lnM' k -^£M' k>ol \nM' k>a +^-^lniV', (101) 



fc=2 



fc=2 



and the complete expression for the annealed entropy is 

sLnlN', {M' k }; {M' k J} = s[[{M' k J] - s' [N', {M' k }}. (102) 
The numerical predictions for the point where the entropy becomes zero are provided in Table ITIl 

VIII. QAA ALGORITHM ON A CORE: EXTENDED LANDSCAPES 

Once we have reexpressed the entropy on a reduced graph, it is only natural to implement the 
quantum adiabatic evolution directly on the reduced graph. Correspondingly, we shall need to 
recompute landscapes for the reduced graph. 

A. Positive K-NAE SAT 

This case entails least difficulty, since we can use same landscape parameters {A4^}, the differ- 
ence being the exclusion of total spin from the list of parameters. The central quantity-£'({.M / cl })- 
is still expressed in a similar form: 



f({KJ) = 5>> kex P 



tr,k 



2 z_, «y dM , a dM ,^ 



(103) 
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Notably, the sum over k is now restricted to |k| ^ 2. The derivatives of the entropy are still 



<*>ann _ ^ ln ^ _ (1Q4) 



dM.' a ^ p 



p 



For a, a' differing in exactly one position, we shall have 



ann ann 



ln^. (105) 



2 V^'„ 9^,7 2 /£' 
Next, using c CT) k = ;§(/•*«) fe< */^c*!> f° r |k| ^ 2 and Eq. d42l . we are able to write, in terms of 
generating function G(fx): 

aiM'j) = \g ( j2 ( l06 ) 

\<a,a'> / 



Substituting M! a f ' jJP a = M' a j \x a , we can rewrite this as 



- gm i G (m-) g yw&ji^ vm ^) ■ <io7) 

Note that for the landscapes on the original graph without spin variable, this expression is valid 
with G{ji) redefined to be G(/i) = e^. 



B. Positive 1-in-if SAT 



The notion of landscape parameters has to be generalized, since clauses of any length can 
appear. Correspondingly, we choose a set {M.' k a } to serve as landscape parameters. In contrast to 
K-NAE SAT, we have a total of £)fL 2 (fc + 1) = (K -1)(K -2)/2 parameters (with symmetries 
of M' ka taken into account). 

The derivation of landscapes can be generalized to include several types of clauses. The final 
answer is itself a generalization of Eq. (l65t 

As before G({fi k }) = exp (j2?=2 ^) ~ e ^ ~ E£s ^ 
We can exploit the symmetry of M.' k a to write 

fc-i 

E \J M K«M> Kct , = E \/ (m + l)i ~ k ' m ^ M UK, m+ r (109) 

<a,a'> m=0 

Same is possible for K-NAE SAT; in fact it is precisely this form that is used in numerical calcu- 
lations. 
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C. Numerical Results 

Here we provide the numerical results for the satisfiability transition as determined by max- 
imizing the entropy for energy E = and solving s^ nn = 0. And we also list the location of 
dynamical transition, indicated by global bifurcation in f — e' — IT. All results are expressed in 
terms of connectivity of the original random graph for easy comparison with Table U 

We observed that this refinement of our analysis leaves j d and 7 C of iT-NAE SAT essentially 
unchanged. This is the manifestation of the fact that if either K or 7 is large, the core is not much 
different from the original graph (i.e. that q 1). In contrast, the difference for 1-in-K SAT 
is quite perceptible. For consistency, we compare the location of dynamic phase transition com- 
puted on the core to that computed on the original graph using only M m as landscape parameters 
(omitting total spin, as it is not among parameters for computations performed on a core). We also 
include new better bounds on static transition. All results are summarized in Table ITTI and Fig. |9] 
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6 


7 


8 
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10 


Improved 


i d 


0.653 


0.609 


0.535 
0.553 


0.469 
0.507 


0.421 
0.468 


0.379 
0.435 


0.344 
0.407 


0.317 
0.382 


Old 


id 




0.671 


0.552 


0.471 


0.413 


0.368 


0.333 


0.304 


7c 


0.805 


0.676 


0.609 


0.548 


0.500 


0.461 


0.428 


0.400 



TABLE II: Dynamic and static transition for positive 1-in-K SAT with improved annealing approximation 
and with old method using {M rn } as landscape parameters. Prediction of 7^ for K = 3 compares favorably 
with result of simulations of j' c 0.63. No value (-) indicates the ebasence of dynamical transition. 

IX. CONCLUSION 

We have formulated an ansatz of landscapes and studied the complexity of the quantum adia- 
batic algorithm within the annealing approximation and found the existence of a dynamic transi- 
tion and a hard(exponential) region above that dynamic transition. However, a similar analysis of 
simulated annealing did not reveal any phase transitions. We explain this as follows. The anneal- 
ing approximation should fail for sufficiently small energies. It is commonly known that simulated 
annealing can find suboptimal solutions with very small energies very efficiently, but it takes an 
exponentially long time to actually reach the ground state. The annealing approximation does not 
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FIG. 9: Static (circles) and dynamic (diamonds) transition for various values of K for positive 1-m-K SAT. 

correctly describe very small energies and cannot be used to establish its complexity. Note that we 
can reconcile this with the fact that the annealing approximation becomes exact in the limit when 
number of bits in a clause K — > oo: if the annealing approximation fails for some E < Ek we 
expect that Ek is decreasing to zero as K increases. However for any finite K, the free energy 
computed within the annealing approximation is free from any singularities indicative of a phase 
transition. To study the complexity of simulated annealing one needs to use the tools of spin glass 



theory, in particular, the replica trick 1 19 



20 



(see also below). 

In contrast, in our analysis of the quantum adiabatic algorithm, we observed a first-order phase 
transition, and, importantly, it happens for energies E* = O(Eoo) (where E^ is the expected 
energy at infinite temperature, E^ = Yl z ^z)- Moreover, the energies on both sides of the 
transition, relative to E^ seem not to change appreciably with increasing K. Since the annealing 
approximation for this range of energies can be used, the prediction for the dynamic transition 
should survive, though the exact numerical values may acquire corrections. We have recomputed 
the dynamic transition with simplified energy-only landscapes (see Fig. fTOl . For 1-in-K SAT one 
can clearly see that the relative correction quickly diminishes. We believe that same happens for K- 
NAE SAT if sufficiently large K's are considered. If this indeed holds, it serves as a corroboration 
that our results are correct numerically for large K. It should be noted that the large-K limit 
corresponds to the so-called random energy model, where one does not expect to perform better 
then 0(2 N ) via any quantum algorithm. 



The idea of using energy-only landscapes was present in as well as Q and |261. A 
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jump in the time-dependence of the expected energy value was seen in numerical simulations 
J3I, indicative of first-order phase transition, though a different ensemble was considered (only 
instances having a unique solution were considered). 

We also attempted to go beyond simple annealing approximation and studied the dynamical 
transition using its refinement. For that we developed a polynomial mapping of the optimization 
problem defined on a full graph onto the problem defined on its subgraph (a core) where disorder- 
related fluctuations are significantly reduced and annealing approximation is expected to perform 
much better. As a test we used annealing approximation on a core to calculate the position 7 C 
of a static (satisfiability) transition where the entropy of the state with E=0 vanishes. We also 
computed 7 C numerically and found it to be very close to the analytical result. We then studied the 
dynamics of quantum adiabatic evolution algorithm on a core using an extended set of landscape 
functions and found that old results obtained on a full graph are reproduced qualitatively. This 
supports our earlier prediction that the location of the phase transition is not very sensitive to exact 
nature of annealing approximation employed. 

0.2 
0.18 
0.16 
0.14 
^ 0.12 

Y 0.1 

±? 0.08 
0.06 
0.04 
0.02 

°456789 10 
K 

FIG. 10: Relative difference between predictions for the dynamical phase transition point in the case of full 
(7d) and energy-only (7^) landscapes vs of K for 1-in-K SAT (circles) and K-NAE SAT (squares). 

We emphasize that different versions of annealing approximation employed in this paper de- 
scribe the phase transition as a global bifurcation between two macroscopic states (pure states) 
in the space of macroscopic variables defined by a set of landscape functions. The complexity 
is due to tunneling between the pure states. In contrast, spin glass theory predicts the existence 
of an infinite number of pure states [23] at sufficiently small energies. On the other hand, as we 




36 



mentioned above the first-order quantum phase transition occurs for large energies E* and this has 
been confirmed with an improved annealing approximation. 

Although the transition is seemingly absent for small K, a better approach (as compared to 
annealing approximation) may reveal it. Moreover, we believe that if this happens the order of the 
transition will remain unchanged, suggesting that the disorder may be irrelevant for the determi- 
nation of the order of the phase transition and, consequently, for the complexity of the quantum 
adiabatic algorithm. That is, the exponential complexity is not due to the true combinatorial com- 
plexity of the underlying random optimization problem but rather due to peculiarities of the driver 
term and a particular ensemble of random instances considered. 

A future extension of the present work is to include sufficiently large (possibly infinite) number 
of landscape parameters, thereby making annealing approximation increasingly precise. In this 
regard we recall that 1 -bit-flip conditional distribution over landscape parameters employed in this 
paper d42b can be expressed via the set of coefficients {c CT k} that are concentrations of binary 
variables in a given string with different types (a, k) of an immediate neighborhood. In fact, 
these coefficients themselves can be used in an extended set of landscape parameters x. Then 
appropriate effective potential /(x, T) can be introduced and its bifurcation can be studied when 
T varies from oo to 0. Furthermore, one can consider introducing progressively larger sets of 
landscape functions by defining neighborhoods of progressively larger size and using the well- 
known property that local structure of random (hyper)graph is tree-like I27I . 
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APPENDIX A: ON THE NP-COMPLETENESS OF POSITIVE 1-IN-K SAT AND 

POSITIVE K-NAE-SAT. 

We set out to prove that both Positive 1-in-K Sat and K-NAE-Sat are NP-complete. It is 
straightforward to see that it takes a polynomial time to verify the assignment, so these problems 
are in NP. We now prove that they are as hard as the Satisfiability problem, which is NP-complete, 
by showing that any boolean formula can be represented as an instance of these. 

1. Positive 1-in-K Sat 

A clause of type (x, . . . , x, y) necessarily implies x = and y — 1; hence we can represent 
constants and I. A clause of type (0, . . . ,0,x,y) implies x = -<y. Finally, a clause of type 
(0, . . . , 0, x, y, z) is equivalent to a 3-clause (x, y, z) so that we can restrict ourselves to K = 3 
without losing generality. 

For K = 3, immediately observe that three clauses (x, z, u')(y, z, u")(u, u', u") with free vari- 
ables u, u', u" implies z = -<(x A y). This basic building block is in fact sufficient to build any 
boolean formula, as a result, any boolean formula can be cast as an 1-in-K SAT formula. 

2. Positive if-NAE-Sat 

A clause of type {x,...,x, y) necessarily implies x = ->y, and (x, . . . ,x, y, z) is equivalent 
to (x,y,z) so we once again restrict ourselves to K = 3. In contrast to l-in-i^ problem, we 
shall require a non-trivial representation of false or true. We will use pairs of variables to denote 
variables of the boolean formula. Pairs 00 or 11 will represent value false and pairs 01 or 10 will 
represent true. 

The next building block, (x, y, t)(y, z, t)(z, x, t) ensures that t — 1 if the majority of x, y, z 
are and t — if the majority are 1. We shall use a shorthand f(t; x, y, z) to denote this. The 
expression f(z\\ xi,yi, y-i)f{z2\ x 2 , 2/1,2/2) then ensures z = x A y where x, y, z are represented 
as pairs x\x 2 , yiy 2 , z\z 2 as indicated above. The operation of negation is trivial to represent: if 
x = x x x 2 then ->x = (~^xi)x 2 . These two are sufficient to construct any boolean formula. 
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APPENDIX B: NEXT ORDER APPROXIMATION FOR LANDSCAPES 

A better approximation for the values of critical clause-to-variable rations can be obtained if 
we specify the constraint that the distribution of vertex degrees be Poisson (as it is supposed to be 
in a random hypergraph To be precise, we specify that 

E c ^ II 6 ( E - k ?) = c iM - rfu e " K - (B1) 

Consequently, with this constraint the following expression for c CT) k is obtained: 

Ca,k = c {k } , (B2) 

Z {k P } 

where we use 

(\ k p / \ k p 

+e- A n + • 

a J p \ a / 

Annealed entropy can be rewritten in the form 

Sann[g, {M^}} = min < -Xq + M a In — ^ + In Z [A, yP a } \ -Y^ M a teM a +'rlnj-jK, 

K P,ct ) a 

(B4) 

where In Z is given by 

ln z = E C <M ln ( B5 ) 

The equations relating q, {M a } and A, are given in (ETt.lKSb. 

Similarly to Sec. IV Bl we will use the notation (!5Ut . Since \nZ depends only on A and /j±, 
dlnZ/d^ depends only on o p . Therefore, M a / ^i p a = M ap /fi ap . Correspondingly, 

s ann [g, {M a }] = min j -\q + M + ln — '- - + M- In — - + ln Z \ 

-^A^a.lnA^« + 7^7-7^ (B6) 

For convenience, we introduce new variables 

^ ± = f xe ±h . (B7) 

We then readily obtain 

ln Z = 7 ln ft + Y c k ln[2 cosh(A + jfe/i)] , (B8) 
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(where k = J2 P K an d Ck = ^fp e 7i ^) and /i drops out of the expression for s ann altogether: 



Sann[g,{A^«}] = min |-Ag - (M + - M-)h + ^c fe ln[2cosh(A + fc/t)]| (B9) 
- y^ j M a \nM tx + M+\nM+ + M-\nM- + 7I117- 7K 



It is easy to see from this expression what the equations for A, /i are: 

c k tanh[A + kh] = q, 

k 

&c fc tanh[A + kh] = M+-M-. 



(BIO) 



We now turn our attention to the function £(q, {M. a }) given by (l42l) with 9 and y a evaluated from 
Eqs. (l56T).d58t. The computation of e Va ~ y <*' yields 



ya-y a / 




(Bll) 



Multiplied by pP a this becomes 



e v<*-y a > 



VMoM< 



r/t. 



y/M+M- 

The expression for £(q, {M a }) can be written in the form (cf. (I^Ol ) 



(B12) 



{^}) = Ey^E eAs+ " s Il 



(B13) 



with the internal sum running over k consistent with a set of {k p } (IB 1 b . Substituting the quantities 
defined above this becomes 

i(q,{M a }) = 5> {M ^ ^ '—. (B14) 

{k P } 



J {k P } 



After some transformations we finally obtain 



{AU) = I> 



cosh[A + kh] 



(B15) 



where A, h are given by (l50l) . (IB10l) . Using symmetric ansatz (l66l) it is straightforward to calculate 
from dB15l> the restricted function £(q, {M m }) (cf. (loTt). We must note however, that although 
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this represents a next-order improvement over annealing approximation, the relative changes in y c 
and 7 d computed with this improved approximation are nearly imperceptible (~ 10~ 4 ). 
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